

************************
** Study 2: Figure S2.16
************************



** 2017 Results
***************

* Covariates macro
global covars_binaries = "gender1 gender2 gender3 agecat1 agecat2 agecat3 agecat4 agecat5 agecat6 edu1 edu2 edu3 edu4 edu5 edu6 pol_interest1 pol_interest2 pol_interest3 pol_interest4 lr0 lr1 lr2 lr3 lr4 lr5 lr6 lr7 lr8 lr9 lr10 vote2015_1 vote2015_2 vote2015_3 vote2015_4 vote2015_5 vote2015_6 vote2015_7 vote2015_8 vote2015_9 vote2015_10 refvote2016_1 refvote2016_2 refvote2016_3 region1 region2 region3"

* Open dataset
use "LondonBridgeAttack2017.dta", replace

* Regressions
foreach dv of varlist sec imm britid englid {
reg `dv' postattack $covars_binaries if statadate == td(02jun2017) | statadate == td(04jun2017)
estimates store `dv'1
reg `dv' postattack $covars_binaries if inrange(statadate,td(01jun2017),td(02jun2017)) | inrange(statadate,td(04jun2017),td(05jun2017))
estimates store `dv'2
reg `dv' postattack $covars_binaries if inrange(statadate,td(31may2017),td(02jun2017)) | inrange(statadate,td(04jun2017),td(06jun2017))
estimates store `dv'3
}

* Results table
gen days = _n in 1/3
foreach dv of varlist sec imm britid englid {
gen n_`dv' = .
gen pe_`dv' = .
gen se_`dv' = .
gen p_`dv' = .
gen ll_`dv' = .
gen ul_`dv' = .
sum `dv'
gen sd_`dv' = `r(sd)'
gen d_`dv' = .
gen dlow_`dv' = .
gen dup_`dv' = .
forvalues i = 1/3 {
estimates restore `dv'`i'
replace n_`dv' = `e(N)' if days == `i'
estimates replay `dv'`i'
matrix R = r(table)
replace pe_`dv' = R[1,1] if days == `i'
replace se_`dv' = R[2,1] if days == `i'
replace p_`dv' = R[4,1] if days == `i'
replace ll_`dv' = R[5,1] if days == `i'
replace ul_`dv' = R[6,1] if days == `i'
replace d_`dv' = pe_`dv'/sd_`dv' if days == `i'
replace dlow_`dv' = pe_`dv'/sd_`dv' - 1.96*se_`dv'/sd_`dv' if days == `i'
replace dup_`dv' = pe_`dv'/sd_`dv' + 1.96*se_`dv'/sd_`dv' if days == `i'
}
}

drop newid - referrer
drop _est*
keep in 1/3
gen year = 2017
order year

save "MainResults2017.dta", replace	




** 2019 Results
***************

* Covariates macro
global covars_binaries = "gender1 gender2 gender3 agecat1 agecat2 agecat3 agecat4 agecat5 agecat6 edu1 edu2 edu3 edu4 edu5 edu6 pol_interest1 pol_interest2 pol_interest3 lr0 lr1 lr2 lr3 lr4 lr5 lr6 lr7 lr8 lr9 lr10 vote2017_1 vote2017_2 vote2017_3 vote2017_4 vote2017_5 vote2017_6 vote2017_7 vote2017_8 vote2017_9 vote2017_10 vote2017_11 refvote2016_1 refvote2016_2 refvote2016_3 refvote2016_4 region1 region2 region3"

* Open dataset
use "LondonBridgeAttack2019.dta", replace

* One-day comparison
********************

* How many observations to add
local toadd_sec = round((6323 - 3513) / 2)
local toadd_imm = round((6274 - 3454) / 2) 
local toadd_britid = round((5908 - 3080) / 2) 
local toadd_englid = round((3880 - 2427) / 2) 

* Identify control observations
gen control = (statadate == td(28nov2019))

foreach dv of varlist sec imm britid englid {

* Actual results
reg `dv' postattack $covars_binaries if statadate == td(28nov2019) | statadate == td(30nov2019)
estimates store `dv'1

* Add artificial control observations
count
local addcontrol = `r(N)' + `toadd_`dv''
set obs `addcontrol'
gen newcontrol = 1 if id == .
foreach var of varlist `dv' postattack $covars_binaries {
	sum `var' if control == 1
	replace `var' = `r(mean)' if newcontrol == 1
}

* Add artificial treatment observations
count
local addtreat = `r(N)' + `toadd_`dv''
set obs `addtreat'
gen newtreat = 1 if id == . & newcontrol == .
replace postattack = 1 if newtreat == 1
foreach var of varlist `dv' $covars_binaries {
	sum `var' if control == 1 & newcontrol != 1
	replace `var' = `r(mean)' if newtreat == 1
}

* Lower bound estimate (zero treatment effect for artificial)
reg `dv' postattack $covars_binaries if statadate == td(28nov2019) | statadate == td(30nov2019) | newcontrol == 1 | newtreat == 1
estimates store low_`dv'1

* Medium low bound estimate (50% treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'1
replace `dv' = `mean' + 0.5 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if statadate == td(28nov2019) | statadate == td(30nov2019) | newcontrol == 1 | newtreat == 1
estimates store medlow_`dv'1

* Medium upper bound estimate (150% treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'1
replace `dv' = `mean' + 1.5 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if statadate == td(28nov2019) | statadate == td(30nov2019) | newcontrol == 1 | newtreat == 1
estimates store medup_`dv'1

* Upper bound estimate (2x treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'1
replace `dv' = `mean' + 2 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if statadate == td(28nov2019) | statadate == td(30nov2019) | newcontrol == 1 | newtreat == 1
estimates store up_`dv'1

* Drop artificial observations
drop if newcontrol == 1 | newtreat == 1
drop newcontrol newtreat
}

* Drop control identifier
drop control




* Two-day comparison
********************

* How many observations to add
local toadd_sec = round((12511 - 8575) / 2)
local toadd_imm = round((12424 - 8481) / 2) 
local toadd_britid = round((11657 - 7483) / 2) 
local toadd_englid = round((7796 - 6038) / 2) 

* Identify control observations
gen control = (statadate == td(28nov2019))
replace control = 1 if statadate == td(27nov2019)

foreach dv of varlist sec imm britid englid {

* Actual results
reg `dv' postattack $covars_binaries if inrange(statadate,td(27nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(1dec2019))
estimates store `dv'2

* Add artificial control observations
count
local addcontrol = `r(N)' + `toadd_`dv''
set obs `addcontrol'
gen newcontrol = 1 if id == .
foreach var of varlist `dv' postattack $covars_binaries {
	sum `var' if control == 1
	replace `var' = `r(mean)' if newcontrol == 1
}

* Add artificial treatment observations
count
local addtreat = `r(N)' + `toadd_`dv''
set obs `addtreat'
gen newtreat = 1 if id == . & newcontrol == .
replace postattack = 1 if newtreat == 1
foreach var of varlist `dv' $covars_binaries {
	sum `var' if control == 1 & newcontrol != 1
	replace `var' = `r(mean)' if newtreat == 1
}

* Lower bound estimate (zero treatment effect for artificial)
reg `dv' postattack $covars_binaries if inrange(statadate,td(27nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(1dec2019)) | newcontrol == 1 | newtreat == 1
estimates store low_`dv'2

* Medium low bound estimate (50% treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'2
replace `dv' = `mean' + 0.5 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if inrange(statadate,td(27nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(1dec2019)) | newcontrol == 1 | newtreat == 1
estimates store medlow_`dv'2

* Medium upper bound estimate (150% treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'2
replace `dv' = `mean' + 1.5 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if inrange(statadate,td(27nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(1dec2019)) | newcontrol == 1 | newtreat == 1
estimates store medup_`dv'2

* Upper bound estimate (2x treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'2
replace `dv' = `mean' + 2 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if inrange(statadate,td(27nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(1dec2019)) | newcontrol == 1 | newtreat == 1
estimates store up_`dv'2

* Drop artificial observations
drop if newcontrol == 1 | newtreat == 1
drop newcontrol newtreat
}

* Drop control identifier
drop control







* Three-day comparison
********************

* How many observations to add
local toadd_sec = round((22240 - 12923) / 2)
local toadd_imm = round((22077 - 12794) / 2) 
local toadd_britid = round((20660 - 11227) / 2) 
local toadd_englid = round((13745 - 9086) / 2) 

* Identify control observations
gen control = (statadate == td(28nov2019))
replace control = 1 if statadate == td(27nov2019)
replace control = 1 if statadate == td(26nov2019)


foreach dv of varlist sec imm britid englid {

* Actual results
reg `dv' postattack $covars_binaries if inrange(statadate,td(26nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(2dec2019))
estimates store `dv'3


* Add artificial control observations
count
local addcontrol = `r(N)' + `toadd_`dv''
set obs `addcontrol'
gen newcontrol = 1 if id == .
foreach var of varlist `dv' postattack $covars_binaries {
	sum `var' if control == 1
	replace `var' = `r(mean)' if newcontrol == 1
}

* Add artificial treatment observations
count
local addtreat = `r(N)' + `toadd_`dv''
set obs `addtreat'
gen newtreat = 1 if id == . & newcontrol == .
replace postattack = 1 if newtreat == 1
foreach var of varlist `dv' $covars_binaries {
	sum `var' if control == 1 & newcontrol != 1
	replace `var' = `r(mean)' if newtreat == 1
}


* Lower bound estimate (zero treatment effect for artificial)
reg `dv' postattack $covars_binaries if inrange(statadate,td(26nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(2dec2019)) | newcontrol == 1 | newtreat == 1
estimates store low_`dv'3

* Medium low bound estimate (50% treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'3
replace `dv' = `mean' + 0.5 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if inrange(statadate,td(26nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(2dec2019)) | newcontrol == 1 | newtreat == 1
estimates store medlow_`dv'3

* Medium upper bound estimate (150% treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'3
replace `dv' = `mean' + 1.5 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if inrange(statadate,td(26nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(2dec2019)) | newcontrol == 1 | newtreat == 1
estimates store medup_`dv'3

* Upper bound estimate (2x treatment effect for artificial)
sum `dv' if control == 1 & newcontrol != 1
local mean = `r(mean)'
estimates restore `dv'3
replace `dv' = `mean' + 2 * _b[postattack] if newtreat == 1
reg `dv' postattack $covars_binaries if inrange(statadate,td(26nov2019),td(28nov2019)) | inrange(statadate,td(30nov2019),td(2dec2019)) | newcontrol == 1 | newtreat == 1
estimates store up_`dv'3

* Drop artificial observations
drop if newcontrol == 1 | newtreat == 1
drop newcontrol newtreat
}

* Drop control identifier
drop control





* Results table
gen estimate = "Actual" in 1/3
replace estimate = "Low" in 4/6
replace estimate = "MedLow" in 7/9
replace estimate = "MedHigh" in 10/12
replace estimate = "High" in 13/15
gen days = _n in 1/3
replace days = _n - 3 in 4/6
replace days = _n - 6 in 7/9
replace days = _n - 9 in 10/12
replace days = _n - 12 in 13/15
foreach dv of varlist sec imm britid englid {
gen n_`dv' = .
gen pe_`dv' = .
gen se_`dv' = .
gen p_`dv' = .
gen ll_`dv' = .
gen ul_`dv' = .
sum `dv'
gen sd_`dv' = `r(sd)'
gen d_`dv' = .
gen dlow_`dv' = .
gen dup_`dv' = .
forvalues i = 1/3 {
	* Actual results
estimates restore `dv'`i'
replace n_`dv' = `e(N)' if days == `i' & estimate == "Actual"
estimates replay `dv'`i'
matrix R = r(table)
replace pe_`dv' = R[1,1] if days == `i' & estimate == "Actual"
replace se_`dv' = R[2,1] if days == `i' & estimate == "Actual"
replace p_`dv' = R[4,1] if days == `i' & estimate == "Actual"
replace ll_`dv' = R[5,1] if days == `i' & estimate == "Actual"
replace ul_`dv' = R[6,1] if days == `i' & estimate == "Actual"
	* Lower bound
estimates restore low_`dv'`i'
replace n_`dv' = `e(N)' if days == `i' & estimate == "Low"
estimates replay low_`dv'`i'
matrix R = r(table)
replace pe_`dv' = R[1,1] if days == `i' & estimate == "Low"
replace se_`dv' = R[2,1] if days == `i' & estimate == "Low"
replace p_`dv' = R[4,1] if days == `i' & estimate == "Low"
replace ll_`dv' = R[5,1] if days == `i' & estimate == "Low"
replace ul_`dv' = R[6,1] if days == `i' & estimate == "Low"
	* Medium lower bound
estimates restore medlow_`dv'`i'
replace n_`dv' = `e(N)' if days == `i' & estimate == "MedLow"
estimates replay medlow_`dv'`i'
matrix R = r(table)
replace pe_`dv' = R[1,1] if days == `i' & estimate == "MedLow"
replace se_`dv' = R[2,1] if days == `i' & estimate == "MedLow"
replace p_`dv' = R[4,1] if days == `i' & estimate == "MedLow"
replace ll_`dv' = R[5,1] if days == `i' & estimate == "MedLow"
replace ul_`dv' = R[6,1] if days == `i' & estimate == "MedLow"

	* Medium upper bound
estimates restore medup_`dv'`i'
replace n_`dv' = `e(N)' if days == `i' & estimate == "MedHigh"
estimates replay medup_`dv'`i'
matrix R = r(table)
replace pe_`dv' = R[1,1] if days == `i' & estimate == "MedHigh"
replace se_`dv' = R[2,1] if days == `i' & estimate == "MedHigh"
replace p_`dv' = R[4,1] if days == `i' & estimate == "MedHigh"
replace ll_`dv' = R[5,1] if days == `i' & estimate == "MedHigh"
replace ul_`dv' = R[6,1] if days == `i' & estimate == "MedHigh"

	* Upper bound
estimates restore up_`dv'`i'
replace n_`dv' = `e(N)' if days == `i' & estimate == "High"
estimates replay up_`dv'`i'
matrix R = r(table)
replace pe_`dv' = R[1,1] if days == `i' & estimate == "High"
replace se_`dv' = R[2,1] if days == `i' & estimate == "High"
replace p_`dv' = R[4,1] if days == `i' & estimate == "High"
replace ll_`dv' = R[5,1] if days == `i' & estimate == "High"
replace ul_`dv' = R[6,1] if days == `i' & estimate == "High"
	* Cohen's d
replace d_`dv' = pe_`dv'/sd_`dv' if days == `i'
replace dlow_`dv' = pe_`dv'/sd_`dv' - 1.96*se_`dv'/sd_`dv' if days == `i'
replace dup_`dv' = pe_`dv'/sd_`dv' + 1.96*se_`dv'/sd_`dv' if days == `i'
}
}

drop id - region3
drop _est*
keep in 1/15
gen year = 2019
order year

save "SensR32019.dta", replace





** Figure
*********


* Bring data in correct shape - 2017
use "MainResults2017.dta", replace
set obs 15
sum year
replace year = `r(mean)'
gen row = _n
gen outcome = ""
gen pe = .
gen ll = .
gen ul = .
gen n = .
forvalues i = 1/3 {
	replace outcome = "sec" if row == `i' + 3
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 3
	sum d_sec if row == `i'
	replace pe = `r(mean)' if row == `i' + 3
	sum dlow_sec if row == `i'
	replace ll = `r(mean)' if row == `i' + 3
	sum dup_sec if row == `i'
	replace ul = `r(mean)' if row == `i' + 3
	sum n_sec if row == `i'
	replace n = `r(mean)' if row == `i' + 3
}
forvalues i = 1/3 {
	replace outcome = "imm" if row == `i' + 6
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 6
	sum d_imm if row == `i'
	replace pe = `r(mean)' if row == `i' + 6
	sum dlow_imm if row == `i'
	replace ll = `r(mean)' if row == `i' + 6
	sum dup_imm if row == `i'
	replace ul = `r(mean)' if row == `i' + 6
	sum n_imm if row == `i'
	replace n = `r(mean)' if row == `i' + 6
}
forvalues i = 1/3 {
	replace outcome = "britid" if row == `i' + 9
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 9
	sum d_britid if row == `i'
	replace pe = `r(mean)' if row == `i' + 9
	sum dlow_britid if row == `i'
	replace ll = `r(mean)' if row == `i' + 9
	sum dup_britid if row == `i'
	replace ul = `r(mean)' if row == `i' + 9
	sum n_britid if row == `i'
	replace n = `r(mean)' if row == `i' + 9
}
forvalues i = 1/3 {
	replace outcome = "englid" if row == `i' + 12
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 12
	sum d_englid if row == `i'
	replace pe = `r(mean)' if row == `i' + 12
	sum dlow_englid if row == `i'
	replace ll = `r(mean)' if row == `i' + 12
	sum dup_englid if row == `i'
	replace ul = `r(mean)' if row == `i' + 12
	sum n_englid if row == `i'
	replace n = `r(mean)' if row == `i' + 12
}
drop in 1/3
drop n_sec - row
save temp2017.dta, replace




* Bring results in correct shape - 2019
use "SensR32019.dta", replace
set obs 15
sum year
replace year = `r(mean)'
gen row = _n
gen outcome = ""
gen pe = .
gen ll = .
gen ul = .
gen n = .
gen low_pe = .
gen low_ll = .
gen low_ul = .
gen low_n = .
gen medlow_pe = .
gen medlow_ll = .
gen medlow_ul = .
gen medlow_n = .
gen medup_pe = .
gen medup_ll = .
gen medup_ul = .
gen medup_n = .
gen up_pe = .
gen up_ll = .
gen up_ul = .
gen up_n = .

forvalues i = 1/3 {
	replace outcome = "sec" if row == `i' + 3
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 3
	
	sum d_sec if row == `i' & estimate == "Actual"
	replace pe = `r(mean)' if row == `i' + 3
	sum dlow_sec if row == `i' & estimate == "Actual"
	replace ll = `r(mean)' if row == `i' + 3
	sum dup_sec if row == `i' & estimate == "Actual"
	replace ul = `r(mean)' if row == `i' + 3
	sum n_sec if row == `i' & estimate == "Actual"
	replace n = `r(mean)' if row == `i' + 3
	
	sum d_sec if row == `i' + 3 & estimate == "Low"
	replace low_pe = `r(mean)' if row == `i' + 3
	sum dlow_sec if row == `i' + 3 & estimate == "Low"
	replace low_ll = `r(mean)' if row == `i' + 3
	sum dup_sec if row == `i' + 3 & estimate == "Low"
	replace low_ul = `r(mean)' if row == `i' + 3
	sum n_sec if row == `i' + 3 & estimate == "Low"
	replace low_n = `r(mean)' if row == `i' + 3
	
	sum d_sec if row == `i' + 6 & estimate == "MedLow"
	replace medlow_pe = `r(mean)' if row == `i' + 3
	sum dlow_sec if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ll = `r(mean)' if row == `i' + 3
	sum dup_sec if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ul = `r(mean)' if row == `i' + 3
	sum n_sec if row == `i' + 6 & estimate == "MedLow"
	replace medlow_n = `r(mean)' if row == `i' + 3	
	
	sum d_sec if row == `i' + 9 & estimate == "MedHigh"
	replace medup_pe = `r(mean)' if row == `i' + 3
	sum dlow_sec if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ll = `r(mean)' if row == `i' + 3
	sum dup_sec if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ul = `r(mean)' if row == `i' + 3
	sum n_sec if row == `i' + 9 & estimate == "MedHigh"
	replace medup_n = `r(mean)' if row == `i' + 3	
	
	sum d_sec if row == `i' + 12 & estimate == "High"
	replace up_pe = `r(mean)' if row == `i' + 3
	sum dlow_sec if row == `i' + 12 & estimate == "High"
	replace up_ll = `r(mean)' if row == `i' + 3
	sum dup_sec if row == `i' + 12 & estimate == "High"
	replace up_ul = `r(mean)' if row == `i' + 3
	sum n_sec if row == `i' + 12 & estimate == "High"
	replace up_n = `r(mean)' if row == `i' + 3		
}

forvalues i = 1/3 {
	replace outcome = "imm" if row == `i' + 6
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 6
		
	sum d_imm if row == `i' & estimate == "Actual"
	replace pe = `r(mean)' if row == `i' + 6
	sum dlow_imm if row == `i' & estimate == "Actual"
	replace ll = `r(mean)' if row == `i' + 6
	sum dup_imm if row == `i' & estimate == "Actual"
	replace ul = `r(mean)' if row == `i' + 6
	sum n_imm if row == `i' & estimate == "Actual"
	replace n = `r(mean)' if row == `i' + 6
	
	sum d_imm if row == `i' + 3 & estimate == "Low"
	replace low_pe = `r(mean)' if row == `i' + 6
	sum dlow_imm if row == `i' + 3 & estimate == "Low"
	replace low_ll = `r(mean)' if row == `i' + 6
	sum dup_imm if row == `i' + 3 & estimate == "Low"
	replace low_ul = `r(mean)' if row == `i' + 6
	sum n_imm if row == `i' + 3 & estimate == "Low"
	replace low_n = `r(mean)' if row == `i' + 6
	
	sum d_imm if row == `i' + 6 & estimate == "MedLow"
	replace medlow_pe = `r(mean)' if row == `i' + 6
	sum dlow_imm if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ll = `r(mean)' if row == `i' + 6
	sum dup_imm if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ul = `r(mean)' if row == `i' + 6
	sum n_imm if row == `i' + 6 & estimate == "MedLow"
	replace medlow_n = `r(mean)' if row == `i' + 6
	
	sum d_imm if row == `i' + 9 & estimate == "MedHigh"
	replace medup_pe = `r(mean)' if row == `i' + 6
	sum dlow_imm if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ll = `r(mean)' if row == `i' + 6
	sum dup_imm if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ul = `r(mean)' if row == `i' + 6
	sum n_imm if row == `i' + 9 & estimate == "MedHigh"
	replace medup_n = `r(mean)' if row == `i' + 6	
	
	sum d_imm if row == `i' + 12 & estimate == "High"
	replace up_pe = `r(mean)' if row == `i' + 6
	sum dlow_imm if row == `i' + 12 & estimate == "High"
	replace up_ll = `r(mean)' if row == `i' + 6
	sum dup_imm if row == `i' + 12 & estimate == "High"
	replace up_ul = `r(mean)' if row == `i' + 6
	sum n_imm if row == `i' + 12 & estimate == "High"
	replace up_n = `r(mean)' if row == `i' + 6			
}

forvalues i = 1/3 {
	replace outcome = "britid" if row == `i' + 9
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 9
	
	sum d_britid if row == `i' & estimate == "Actual"
	replace pe = `r(mean)' if row == `i' + 9
	sum dlow_britid if row == `i' & estimate == "Actual"
	replace ll = `r(mean)' if row == `i' + 9
	sum dup_britid if row == `i' & estimate == "Actual"
	replace ul = `r(mean)' if row == `i' + 9
	sum n_britid if row == `i' & estimate == "Actual"
	replace n = `r(mean)' if row == `i' + 9	
	
	sum d_britid if row == `i' + 3 & estimate == "Low"
	replace low_pe = `r(mean)' if row == `i' + 9
	sum dlow_britid if row == `i' + 3 & estimate == "Low"
	replace low_ll = `r(mean)' if row == `i' + 9
	sum dup_britid if row == `i' + 3 & estimate == "Low"
	replace low_ul = `r(mean)' if row == `i' + 9
	sum n_britid if row == `i' + 3 & estimate == "Low"
	replace low_n = `r(mean)' if row == `i' + 9
	
	sum d_britid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_pe = `r(mean)' if row == `i' + 9
	sum dlow_britid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ll = `r(mean)' if row == `i' + 9
	sum dup_britid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ul = `r(mean)' if row == `i' + 9
	sum n_britid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_n = `r(mean)' if row == `i' + 9	
	
	sum d_britid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_pe = `r(mean)' if row == `i' + 9
	sum dlow_britid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ll = `r(mean)' if row == `i' + 9
	sum dup_britid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ul = `r(mean)' if row == `i' + 9
	sum n_britid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_n = `r(mean)' if row == `i' + 9	
	
	sum d_britid if row == `i' + 12 & estimate == "High"
	replace up_pe = `r(mean)' if row == `i' + 9
	sum dlow_britid if row == `i' + 12 & estimate == "High"
	replace up_ll = `r(mean)' if row == `i' + 9
	sum dup_britid if row == `i' + 12 & estimate == "High"
	replace up_ul = `r(mean)' if row == `i' + 9
	sum n_britid if row == `i' + 12 & estimate == "High"
	replace up_n = `r(mean)' if row == `i' + 9
}

forvalues i = 1/3 {
	replace outcome = "englid" if row == `i' + 12
	sum days if row == `i'
	replace days = `r(mean)' if row == `i' + 12
	
	sum d_englid if row == `i' & estimate == "Actual"
	replace pe = `r(mean)' if row == `i' + 12
	sum dlow_englid if row == `i' & estimate == "Actual"
	replace ll = `r(mean)' if row == `i' + 12
	sum dup_englid if row == `i' & estimate == "Actual"
	replace ul = `r(mean)' if row == `i' + 12
	sum n_englid if row == `i' & estimate == "Actual"
	replace n = `r(mean)' if row == `i' + 12	
	
	sum d_englid if row == `i' + 3 & estimate == "Low"
	replace low_pe = `r(mean)' if row == `i' + 12
	sum dlow_englid if row == `i' + 3 & estimate == "Low"
	replace low_ll = `r(mean)' if row == `i' + 12
	sum dup_englid if row == `i' + 3 & estimate == "Low"
	replace low_ul = `r(mean)' if row == `i' + 12
	sum n_englid if row == `i' + 3 & estimate == "Low"
	replace low_n = `r(mean)' if row == `i' + 12
	
	sum d_englid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_pe = `r(mean)' if row == `i' + 12
	sum dlow_englid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ll = `r(mean)' if row == `i' + 12
	sum dup_englid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_ul = `r(mean)' if row == `i' + 12
	sum n_englid if row == `i' + 6 & estimate == "MedLow"
	replace medlow_n = `r(mean)' if row == `i' + 12
	
	sum d_englid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_pe = `r(mean)' if row == `i' + 12
	sum dlow_englid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ll = `r(mean)' if row == `i' + 12
	sum dup_englid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_ul = `r(mean)' if row == `i' + 12
	sum n_englid if row == `i' + 9 & estimate == "MedHigh"
	replace medup_n = `r(mean)' if row == `i' + 12	
	
	sum d_englid if row == `i' + 12 & estimate == "High"
	replace up_pe = `r(mean)' if row == `i' + 12
	sum dlow_englid if row == `i' + 12 & estimate == "High"
	replace up_ll = `r(mean)' if row == `i' + 12
	sum dup_englid if row == `i' + 12 & estimate == "High"
	replace up_ul = `r(mean)' if row == `i' + 12
	sum n_englid if row == `i' + 12 & estimate == "High"
	replace up_n = `r(mean)' if row == `i' + 12		
}
drop in 1/3
drop estimate n_sec - row
save temp2019.dta, replace



* Merge 2017 and 2019
use temp2017.dta, replace
append using temp2019.dta

* Graph preparations
gen ntext = "{it:N} = "
tostring n, replace
replace ntext = ntext + n if year == 2017
tostring low_n, replace
replace ntext = ntext + low_n if year == 2019
gen ypos = .
replace ypos = 0.70 if days == 1
replace ypos = 0.50 if days == 2
replace ypos = 0.30 if days == 3
replace ypos = ypos + 3 if outcome == "sec"
replace ypos = ypos + 2 if outcome == "imm"
replace ypos = ypos + 1 if outcome == "britid"
gen ypos_actual = ypos
gen ypos_low = ypos + 0.08 if year == 2019
gen ypos_medlow = ypos + 0.04 if year == 2019
gen ypos_medup = ypos - 0.04 if year == 2019
gen ypos_up = ypos - 0.08 if year == 2019
gen x = 0.30


* Graph
twoway ///
	(scatter ypos_low low_pe if days == 1, mcolor(midblue) msymbol(o) msize(medium)) ///
	(rspike low_ll low_ul ypos_low if days == 1, lcolor(midblue) horizontal) ///
	(scatter ypos_low low_pe if days == 2, mcolor(midblue) msymbol(s) msize(medium)) ///
	(rspike low_ll low_ul ypos_low if days == 2, lcolor(midblue) horizontal) ///
	(scatter ypos_low low_pe if days == 3, mcolor(midblue) msymbol(t) msize(medium)) ///
	(rspike low_ll low_ul ypos_low if days == 3, lcolor(midblue) horizontal) ///	
	 ///
	(scatter ypos_medlow medlow_pe if days == 1, mcolor(midblue*.5) msymbol(o) msize(medium)) ///
	(rspike medlow_ll medlow_ul ypos_medlow if days == 1, lcolor(midblue*.5) horizontal) ///
	(scatter ypos_medlow medlow_pe if days == 2, mcolor(midblue*.5) msymbol(s) msize(medium)) ///
	(rspike medlow_ll medlow_ul ypos_medlow if days == 2, lcolor(midblue*.5) horizontal) ///
	(scatter ypos_medlow medlow_pe if days == 3, mcolor(midblue*.5) msymbol(t) msize(medium)) ///
	(rspike medlow_ll medlow_ul ypos_medlow if days == 3, lcolor(midblue*.5) horizontal) ///
	 ///	
	(scatter ypos_medup medup_pe if days == 1, mcolor(red*.5) msymbol(o) msize(medium)) ///
	(rspike medup_ll medup_ul ypos_medup if days == 1, lcolor(red*.5) horizontal) ///
	(scatter ypos_medup medup_pe if days == 2, mcolor(red*.5) msymbol(s) msize(medium)) ///
	(rspike medup_ll medup_ul ypos_medup if days == 2, lcolor(red*.5) horizontal) ///
	(scatter ypos_medup medup_pe if days == 3, mcolor(red*.5) msymbol(t) msize(medium)) ///
	(rspike medup_ll medup_ul ypos_medup if days == 3, lcolor(red*.5) horizontal) ///
	 ///	 
	(scatter ypos_up up_pe if days == 1, mcolor(red*.9) msymbol(o) msize(medium)) ///
	(rspike up_ll up_ul ypos_up if days == 1, lcolor(red*.9) horizontal) ///
	(scatter ypos_up up_pe if days == 2, mcolor(red*.9) msymbol(s) msize(medium)) ///
	(rspike up_ll up_ul ypos_up if days == 2, lcolor(red*.9) horizontal) ///
	(scatter ypos_up up_pe if days == 3, mcolor(red*.9) msymbol(t) msize(medium)) ///
	(rspike up_ll up_ul ypos_up if days == 3, lcolor(red*.9) horizontal) ///
	 ///
	(scatter ypos_actual pe if days == 1, mcolor(gs12) msymbol(o) msize(medium)) ///
	(rspike ll ul ypos_actual if days == 1, lcolor(gs12) horizontal) ///
	(scatter ypos_actual pe if days == 2, mcolor(gs7) msymbol(s) msize(medium)) ///
	(rspike ll ul ypos_actual if days == 2, lcolor(gs7) horizontal) ///
	(scatter ypos_actual pe if days == 3, mcolor(gs2) msymbol(t) msize(medium)) ///
	(rspike ll ul ypos_actual if days == 3, lcolor(gs2) horizontal) ///
	 ///
	(scatter ypos x, msymbol(none) mlabel(ntext) mlabsize(vsmall) mlabangle(horizontal) mlabposition(3)) ///
	, ///
	by(year, noixlabel ixtitle graphregion(fcolor(white) lcolor(white)) bgcolor(white) note("") legend(pos(6))) ///
	ytitle("") yscale(noline range(0 4)) ///
	ylabel(3.5 "Tough security" 2.5 "Anti-immigration" 1.5 "British identity" 0.5 "English identity", angle(horizontal) nogrid) ///
	xtitle("Effect on political attitudes (Cohen's {it:d})", margin(small))  ///
	xline(0, lwidth(thin) lpattern(solid) lcolor(black) extend) ///
	xlabel(-0.10(.10).40,) xmlabel(-0.10(0.05)0.40, ) xscale(noline) ///
	legend(order(25 27 29) label(25 "± 1 day") label(27 "± 2 days") label(29 "± 3 days") rows(1) size(small) keygap(*1) region(lstyle(none) lcolor(white))) ///
	subtitle(, size(large) align(middle) margin(bottom) nobox fcolor(white))  ///
	graphregion(fcolor(white) ifcolor(white) lcolor(white)) plotregion(fcolor(white) lcolor(black)) bgcolor(white) ///
	scheme(s2mono) xsize(12) ysize(7)
gr_edit .b1title.draw_view.setstyle, style(no)
gr_edit .plotregion1.xaxis1[1].style.editstyle majorstyle(tickstyle(show_labels(yes))) editcopy
gr_edit .legend.DragBy 0 12
gr_edit .gmetric_mult = 1.15





** Tidy Up
**********

erase "MainResults2017.dta"
erase "SensR32019.dta"
erase temp2017.dta
erase temp2019.dta

